%Martin M. Andreasen, Jan 2023
addpath(genpath(pwd))
clear
close all
%clc
%% Settings
numCPUs         = feature('numcores');

lambdaSMM        = 1       %Weight in objective function to Shrinkage of based on Method of Moments (SMM)
SMMwithCSreg     = 1       %2 for incl. CS and risk-adj CS moments, 
                           %1 for incl. CS moments,
                           %0 for not including CS moments
inclSurveys      = 0
paramsVersion    = 2       %1 for standard EZ, 
 		                   %2 for model with constant
                           %3 for model with ALFA = 0
                           %4 for model with sigmad = 0
			               %5 for model with constant under Perfect foresight

GH_points       = 5;       %Number of points for numerical integration
xPointsDim      = 0;       %Number of points per state dimension for Accuracy grid. 
                           % If 0, then use sim path
pruning         = 0;       %Compute accuracy without using pruning
mpON            = 1;
orderAppMax     = 5-2;     %Order of approximation for the perturbation solution
nyMacro         = 9;       % Number of "macro equations" - most appear first in y
nyBonds         = 40;	   % Number of "bond equations" most appear last in y

%% Model parameters
% Load information about the model
if numCPUs <= 8
    load([pwd,'\ModelSpecification\setupModel'],'setupModel')
else
    load([pwd,'/ModelSpecification/setupModel'],'setupModel')
end
 
if paramsVersion == 1
    modelName = ['AccuracyGrid0_endYear2019_lambdaSMM',num2str(lambdaSMM),'_SMMwithCSreg',num2str(SMMwithCSreg),'_inclSurveys',num2str(inclSurveys),'_withSE'];
elseif paramsVersion == 2
    modelName = ['AccuracyGrid0_u_endYear2019_lambdaSMM',num2str(lambdaSMM),'_SMMwithCSreg',num2str(SMMwithCSreg),'_inclSurveys',num2str(inclSurveys),'_withSE'];
elseif paramsVersion == 3
    modelName = ['AccuracyGrid0_ALFA0_endYear2019_lambdaSMM',num2str(lambdaSMM),'_SMMwithCSreg',num2str(SMMwithCSreg),'_inclSurveys',num2str(inclSurveys)'];
elseif paramsVersion == 4
    modelName = ['AccuracyGrid_u_SIGMAd0_endYear2019_lambdaSMM1_SMMwithCSreg1_inclSurveys0'];
elseif paramsVersion == 5
    modelName = ['AccuracyGrid_u_PFsolv_endYear2019_lambdaSMM1_SMMwithCSreg1_inclSurveys0'];
    GH_points = 1;  %for perfect foresight solution
end
tmpProj   = load([modelName,'.mat']);
results   = tmpProj.results;
params    = results.model.params;
modelProj2nd = tmpProj.modelYield;

%params.SIGMAmuz = 0;
%results.xGridAccuracy(:,2) = 0;

% Solving model for simulating the economy for this model
if paramsVersion == 5
    % no shocks in perturbation solution 
    params.SIGMAmuz = 0;
    params.SIGMAd = 0;
    params.SIGMAn = 0;
    params.SIGMAp = 0;
    params.SIGMAa = 0;
end
MexOn = 0;
tic
[model,errorMes] = perturbationDSGEmodel(params,orderAppMax,setupModel,MexOn);
toc
   
% Standard deviations of the states in linearized model
nx   = size(model.eta,1);
varX = reshape((eye(nx^2)-kron(model.hx,model.hx))\reshape(model.eta*model.eta',nx^2,1),nx,nx);
stdX = sqrt(max(diag(varX),1D-6));
    
% Setup structure for accuracy
[setupEE.GH_e,setupEE.GH_w] = GaussHermite(GH_points);
setupEE.lowerx     = -3*stdX.*ones(length(model.nx),1);
setupEE.upperx     =  3*stdX.*ones(length(model.nx),1);
setupEE.xPointsDim = xPointsDim;
%setupEE.xSim       = results.xGridAccuracy(1:2000,:)';
setupEE.xSim       = results.xGridAccuracy(:,:)';
setupEE.numCPUs    = numCPUs;

%% Accuracy
extendedPerOn = 0;
RMSE_total = zeros(model.ny,orderAppMax+2);
RMSE_macro = zeros(nyMacro,orderAppMax+2);
RMSE_bonds = zeros(nyBonds,orderAppMax+2);

MAE_total  = zeros(model.ny,orderAppMax+2);
MAE_macro  = zeros(nyMacro,orderAppMax+2);
MAE_bonds  = zeros(nyBonds,orderAppMax+2);

MaxAE_total = zeros(model.ny,orderAppMax+2);
MaxAE_macro = zeros(nyMacro,orderAppMax+2);
MaxAE_bonds = zeros(nyBonds,orderAppMax+2);

meanYields  = zeros(nyBonds,orderAppMax);
stdYields   = zeros(nyBonds,orderAppMax);
for ii= 1:orderAppMax
    ii
    if pruning == 0
        if mpON == 1
            [ResEuler_Per,yGrid,xGrid] = EulerEqError_mp(model,extendedPerOn,ii,setupEE);
        else
            [ResEuler_Per,yGrid,xGrid] = EulerEqError(model,extendedPerOn,ii,setupEE);
        end
    elseif pruning == 1
        if mpON == 1
            [ResEuler_Per,yGrid] = EulerEqErrorPruning_mp(model,ii,setupEE);
        else
            [ResEuler_Per,yGrid] = EulerEqErrorPruning(model,ii,setupEE);
        end
    end
    % Measures of accuracy
    RMSE_total(:,ii)  = sqrt(mean(ResEuler_Per(:,1:nyMacro+nyBonds).^2,1))';
    RMSE_macro(:,ii)  = sqrt(mean(ResEuler_Per(:,1:nyMacro).^2,1))';
    RMSE_bonds(:,ii)  = sqrt(mean(ResEuler_Per(:,nyMacro+1:nyMacro+nyBonds).^2,1))';
    
    MAE_total(:,ii)  = mean(abs(ResEuler_Per(:,1:nyMacro+nyBonds)),1)';
    MAE_macro(:,ii)  = mean(abs(ResEuler_Per(:,1:nyMacro)),1)';
    MAE_bonds(:,ii)  = mean(abs(ResEuler_Per(:,nyMacro+1:nyMacro+nyBonds)),1)';
    
    MaxAE_total(:,ii)  = max(abs(ResEuler_Per(:,1:nyMacro+nyBonds)),[],1)';
    MaxAE_macro(:,ii)  = max(abs(ResEuler_Per(:,1:nyMacro)),[],1)';
    MaxAE_bonds(:,ii)  = max(abs(ResEuler_Per(:,nyMacro+1:nyMacro+nyBonds)),[],1)';
    
    rng(1)
    numSim = 1D5;
    shocks = randn(model.ne,numSim); %Simulating from a random normal
    extendedPerOn      = 0;
    x0                 = model.h0;    % initial state is set at the steady state
    pruningSim         = 1;
    if pruningSim == 0
        [ySim_Per,xSim_Per,xPruning]= simul(x0,shocks,model.g0,model.h0,model.eta,model.derivs,ii,pruningSim);
    elseif pruningSim == 1
        [ySim_Per,xSim_Per,xPruning]= simul(x0,shocks,model.g0,model.h0,model.eta,model.derivs,ii,pruningSim);
    end
    meanYields(:,ii) = mean((-4./(1:40)').*ySim_Per(nyMacro+1:nyMacro+nyBonds,:),2);
    stdYields(:,ii)  = std((-4./(1:40)').*ySim_Per(nyMacro+1:nyMacro+nyBonds,:),0,2);
end

%% Projection at 2nd order
modelProj2ndUsed = modelProj2nd;
% required order for variables: 
% c_cu r_cu pai_cu w_cu l_cu output_cu mc_cu evf_cu vf_cu
% Add short rate to y
modelProj2ndUsed.labely = [modelProj2nd.labely '$r_t$'];
modelProj2ndUsed.gSS    = [modelProj2nd.gSS;modelProj2nd.bySS(1,1)];
modelProj2ndUsed.g0     = [modelProj2nd.g0;modelProj2nd.by0(1,1)];
modelProj2ndUsed.gx     = [modelProj2nd.gx;modelProj2nd.byx(1,:)];
modelProj2ndUsed.Gxx    = [modelProj2nd.Gxx;modelProj2nd.Byxx(1,:)];
modelProj2ndUsed.Gxxx   = [modelProj2nd.Gxxx;modelProj2nd.Byxxx(1,:)];
% Re-order the variables to fit these perturbation codes.
order = zeros(1,nyMacro);
for i=1:nyMacro
    if strcmp(model.labely(1,i),'$mvf_t$')
        tmp = find(strcmp('$vf_t$',modelProj2ndUsed.labely));
    else
        tmp = find(strcmp(model.labely(1,i),modelProj2ndUsed.labely));
    end
    order(1,i) = tmp(1,1);
end
modelProj2ndUsed.labely = modelProj2ndUsed.labely(order);
modelProj2ndUsed.gSS    = modelProj2ndUsed.gSS(order);
modelProj2ndUsed.g0     = modelProj2ndUsed.g0(order);
modelProj2ndUsed.gx     = modelProj2ndUsed.gx(order,:);
modelProj2ndUsed.Gxx    = modelProj2ndUsed.Gxx(order,:);
modelProj2ndUsed.Gxxx   = modelProj2ndUsed.Gxxx(order,:);
modelProj2ndUsed.ny     = length(order)+40;
modelProj2ndUsed.labely
if mpON == 1
    [ResEuler_Proj2nd,xGrid] = EulerEqError_Proj_mp(modelProj2ndUsed,setupEE,2);
else
    [ResEuler_Proj2nd,xGrid] = EulerEqError_Proj(modelProj2ndUsed,setupEE,2);
end
ii = orderAppMax+2;
RMSE_total(:,ii)  = sqrt(mean(ResEuler_Proj2nd(:,1:nyMacro+nyBonds).^2,1))';
RMSE_macro(:,ii)  = sqrt(mean(ResEuler_Proj2nd(:,1:nyMacro).^2,1))';
RMSE_bonds(:,ii)  = sqrt(mean(ResEuler_Proj2nd(:,nyMacro+1:nyMacro+nyBonds).^2,1))';

MAE_total(:,ii)  = mean(abs(ResEuler_Proj2nd(:,1:nyMacro+nyBonds)),1)';
MAE_macro(:,ii)  = mean(abs(ResEuler_Proj2nd(:,1:nyMacro)),1)';
MAE_bonds(:,ii)  = mean(abs(ResEuler_Proj2nd(:,nyMacro+1:nyMacro+nyBonds)),1)';

MaxAE_total(:,ii)  = max(abs(ResEuler_Proj2nd(:,1:nyMacro+nyBonds)),[],1)';
MaxAE_macro(:,ii)  = max(abs(ResEuler_Proj2nd(:,1:nyMacro)),[],1)';
MaxAE_bonds(:,ii)  = max(abs(ResEuler_Proj2nd(:,nyMacro+1:nyMacro+nyBonds)),[],1)';


% Do table            
AccuracyTable = [log10(mean(MAE_macro,1));
                 log10(mean(RMSE_macro,1));
                 log10(mean(MaxAE_macro,1));
                 NaN(1,orderAppMax+2);
                 log10(mean(MAE_bonds,1));
                 log10(mean(RMSE_bonds,1));
                 log10(mean(MaxAE_bonds,1))]    
proj2ndBondsCheck = log10(mean(tmpProj.results.ProjBonds,1))'

modelNameOut = ['Per_Prun',num2str(pruning),'_',modelName];
save([modelNameOut,'.mat'],...
    'RMSE_total','RMSE_macro','RMSE_bonds',...
    'MAE_total','MAE_macro','MAE_bonds',...
    'MaxAE_total','MaxAE_macro','MaxAE_bonds','meanYields','stdYields','params','AccuracyTable');